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Abstract. We introduce a model for temporally disordered directed percolation 
in which the probability of spreading from a vertex (t, x) , where t is the time and 
x is the spatial coordinate, is independent of x but depends on t. Using a very 
efficient algorithm we calculate low-density series for bond percolation on the directed 
square lattice. Analysis of the series yields estimates for the critical point p c and 
various critical exponents which are consistent with a continuous change of the critical 
parameters as the strength of the disorder is increased. 
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1. Introduction 

Directed percolation (DP) can be thought of simply as a percolation process on a directed 
lattice in which connections are allowed only in a preferred direction, e.g., on hyper-cubic 
lattices connections are only allowed along edges connecting vertices with increasing 
coordinates. DP is most commonly interpreted as a growth model and the preferred 
direction t is time. Bond percolation is also a special case of a (d + l)-dimensional 
stochastic cellular automaton On the square lattice the evolution is governed by 
the transition probabilities W(a x \o~i, cr r ), with <Xj = 1 if the vertex i is occupied and 
otherwise, which is the probability of finding the vertex x in state o~ x at time t given 
that, at time t — 1, the vertices x — 1 and x + 1 were in states o~\ and o~ r , respectively. 
Bond percolation corresponds to the choice W(0\ai,a r ) = (1 — p) ai+rTr . The behavior of 
the model is controlled by the spreading probability (or density of bonds) p. When p is 
smaller than a critical value p c all clusters remain finite, and in this sub-critical region 
any initial population will eventually die out. Above p c there is a non-zero probability 
of finding an infinite cluster, the average cluster-size S(p) diverges as p — > p c , and the 
initial population will increase exponentially fast. 

Directed percolation, and a closely related continuous time version the contact 
process jU], serve as the prime examples of models for population growth exhibiting a 
non-equilibrium phase transition into an absorbing state (a state from which the system 
cannot escape), in these cases the state totally devoid of occupied sites. DP type 
transitions are also encountered in many other situations, perhaps most prominently 
in models for chemical reactions [2B E] including heterogeneous catalysis and surface 
reactions [22]- For a recent comprehensive review see ^U] 

Most studies of non-equilibrium systems have been limited to cases without 
disorder, that is the spreading probability is homogeneous in both time and space. 
Obviously in many real systems this idealisation is unrealistic. Often some degree of 
disorder is present, e.g., in a real catalyst some sites may be blocked by impurities 
leading one to consider models with quenched spatial disorder. Likewise in trying to 
model real growth processes one would like to consider temporal disorder so as to model 
changes in conditions from year to year, seasonal changes etc. 

Non-equilibrium models with quenched spatial disorder were considered by Kinzel 
[T7] . He argued that the Harris criterion [8J, well-known from disordered spin-systems, 
should apply for DP as well. Harris argued that in systems with quenched disorder one 
must have dv > 2, where v is the correlation length exponent. This means that quenched 
disorder is a relevant perturbation and thus should change the critical exponents if in 
the pure system dv < 2. The Harris criterion has been rigorously established for a 
large class of disordered systems by Chayes et al [2]. In the case of DP the Harris 
criterion becomes di>± < 2, and quenched disorder is a relevant perturbation for d < 3. 
The effect of quenched disorder was studied numerically by Noest 120] who found 
a marked changed in the static critical exponents independent of the strength of the 
disorder. Moreira and Dickman [T%1 |3] used Monte Carlo simulations to study the site 
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diluted two-dimensional contact process on the square lattice. They found a marked 
change in the static critical exponents j3 and v±_ with a value of u± = 1.00(9) consistent 
with the Harris criterion. The dynamic behaviour was found to be incompatible with the 
usual scaling observed in pure models. In particular critical spreading was found to be 
logarithmic rather than power law, and the survival probability in the sub-critical regime 
decayed algebraically rather that exponentially. This means that the dynamical critical 
correlation length exponent vu is undefined. Janssen investigated the problem 
using renormalisation group theory and confirmed the findings of Moreira and Dickman. 
Recently, Hooyberghs, Igloi and Vanderzande [TTJ [TJj have investigated these problems 
using a real-space renormalisation group framework and the numerical technique of the 
density matrix renormalisation group. They found that for strong enough disorder the 
behaviour is controlled by a strong disorder fixed point, with logarithmic dynamical 
correlations and critical exponents different from the pure system. For weaker disorder 
the numerical evidence was consistent with continuously varying static exponents. 

Kinzel P2j also briefly considered the effect of temporal disorder and argued that a 
Harris type criterion applied but with dv± replaced with un. Since vu = 1.733847(6) JE] 
for (l+l)-dimensional DP this would make temporal disorder a relevant perturbation. 
In a previous paper [T^] we studied temporally disordered directed percolation on the 
square lattice using a simple model in which spreading from some rows was deterministic 
(that is spreading takes place with probability 1) while spreading from the other rows 
took place with probability p as in the pure model. Disorder was introduced by letting 
any given row be "deterministic" with probability a independent of other rows. The 
model was studied using high-density series expansions for the percolation probability 
and Monte Carlo simulations of growing clusters from a single seed. The major finding 
was that the critical exponents changed continuously with the strength of the disorder. 
In particular we found that uu < 2 over a wide range of values in apparent violation of 
the Harris criterion. 

The model used in [T^] is not suitable for a low-density expansion so in this paper 
we study another model of DP with temporal disorder. The spreading probability p(t) 
is chosen at random to be either p or ap with probability |. By varying the parameter a 
we can vary the strength of the disorder. There is an obvious symmetry between the two 
regions a < 1 and a > 1 and we shall therefore only consider the case a > 1. Analysis 
of the series indicate that for any given value of a there is critical point p c {ct) where 
the system has a phase transition. For values of the spreading probability p close to 
p c (oi) the model changes from periods of sub-critical to periods of super-critical growth. 
For large values of a there is a very pronounced difference between the two regimes and 
the population dynamics changes from feast to famine. In the infinite lattice limit the 
system will have arbitrarily long stretches of 'famine' conditions in which the population 
is rapidly reduced and very likely do die out. Note that we can be sure that the model 
will have sub- and super-critical regions. We can simply choose p > p c , where p c is the 
critical point of the pure model, and we will always have a super-critical growth process, 
while for ap < p c the growth process is always sub-critical. 
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In Section |2] we give a description of the method used to derive the low-density 
series for this model. The series is then analysed (for various values of a) and the results 
presented in Section El Finally Section 0] contains a brief summary and discussion of 
the main results. 

2. Calculation of low-density series 

In the low-density phase (p < p c ) many quantities of interest can be derived from the 
pair-connectedness Ct >x {p), which is the probability that the vertex at position x is 
occupied at time t given that the origin was occupied at t = 0. Of particular interest 
are moments of the pair-connectedness 

AVnb) =SZ>* n CWp) (1) 

t X 

Due to symmetry moments involving odd powers of x are identically zero. The remaining 
moments diverge as p approaches p c from below 

/Vn(?>) OC ( Pc - p ^ p - (2) 

In a previous paper [THj we gave a detailed description of how the graph theoretical 
properties of the pair-connectedness PP can be turned into a very efficient algorithm for 
the calculation of low-density series expansions for ordinary directed percolation. The 
series expansions for the disordered model is a simple generalisation of this work. Before 
describing the algorithm for the disorder system we will briefly review the pure case. 

It has been shown pQ that the pair-connectedness can be expressed as a sum over 
all graphs (or finite clusters) formed by taking unions of directed paths connecting the 
origin to the vertex (x,t), 

a>) = 5>(0)pi»i, (3) 

9 

where \g\ is the number of bonds in g. The weight d(g) = (— l) c( - g \ where c(g) is the 
cyclomatic number of the graph g. Note that c(g) is increased by 1 whenever two paths 
join, e.g., if there are two incoming bonds on a vertex. The restriction to unions of 
paths is very strong and one immediate consequence is that graphs with dangling parts 
make no contribution to C^ x an d any contributing graph terminates exactly at (t,x). 
Another way of stating the restriction is that any vertex with an incoming bond must 
have an outgoing bond unless it is the terminal vertex (t,x). Any directed path to a 
vertex whose parallel distance from the origin is t contains at least t bonds. One can 
do much better by using a so-called non-nodal graph expansion [5. to extend the series 
to order 2t. A graph g is nodal if it has a vertex (other than the terminal vertex) 
through which all paths pass. It is clear that each such nodal point effectively works 
as a new origin for the cluster growth, and we can obviously obtain any contributing 
graph by concatenating non-nodal graphs. Note that the graph consisting of a single 
bond is non-nodal so all linear graphs can be obtained by repeated concatenations. 
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Figure 1. Example of how the boundary line (cutting through the filled circles) is 
moved in order to insert another vertex (shaded circle). 

This is the essential idea behind the non-nodal graph expansion, which proceeds in two 
principal steps. First we calculate the contribution C£ x of non-nodal graphs to the 
pair-connectedness. The calculation is done up to a preset order N, where N is limited 
by the available computational resources (primarily physical memory). Next we use 
repeated concatenation operations of C* x to calculate the pair-connectedness C t)X and 
from this we finally calculate various moments p, m , n {p)- 

The calculation of C* x is done efficiently using transfer-matrix techniques. This 
involves drawing a boundary line across a finite slice of the lattice and then moving the 
boundary line such that one adds row after row with each row built up one vertex at a 
time, as illustrated in figure HJ The sum over all contributing graphs is calculated as the 
lattice is constructed. At any given stage the boundary line cuts through a number of, 
say j, vertices. There are two possible states (0 or 1) per vertex, corresponding to vertices 
with (1) or without (0) incoming bonds, leading to a total of 2 J possible boundary 
configurations. However, as explained in not all the 2 J possible configurations are 
required. Many of them can be discarded because they either don't make a contribution 
to C* x or only contributes above order N. The weight of each configuration is given by 
a polynomial P in p truncated at order N. The updating of P, as the boundary line is 
moved to insert a new vertex, depends only on the states of the vertices marked x and 
x + 1 in figure E] and are given simply by 

P(S 1A ) =p 2 P(S lfi )+ P P(S 1A )-p 2 P(S 1A ), 
P(S ,x) = pP(Si, ) + P(S ,i) - pP(Si,x), 
P(S h0 )=pP(S lj0 ), 
P(S ,o)=P(S 0> o). 

Here S a ^ is the configuration before the move which has the vertex at x in state a and 
the vertex at x + 1 in state b, while S a ,b is the configuration after the move which has 
the new vertex at x' in state a and the vertex at x + 1 in state b. For a derivation of 
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the updating rules see [TB] . 

Limiting the calculation to non-nodal contributions is very simple; whenever the 
boundary line reaches the horizontal position one just sets to zero the polynomials of 
states with a single incoming bond. This obviously ensures that configurations with a 
nodal point are deleted from the calculation. The pair-connectedness at the following 
time can be calculated from the states with incoming bonds at nearest neighbour sites 
and no incoming bonds on any other sites. In this way we calculate the two non-nodal 
series C?(p) = E.C*» and X* t {p) = 

The generalisation to DP with temporal disorder is quite simple. First of all the 
configuration weights P will depend on two variables u and v corresponding to the 
two spreading probabilities. Lets consider the situation in which a row has just been 
completed. Since the probability is | that spreading from a given row occurs with 
probability u or v, respectively, obviously P(u, v) = P(v, u). Next we insert a new row. 
Formally we can express the weight of any new configuration S as a weighted sum over 
the weights of the configurations S in the row above 

Ps(u, v) = W ( u ) p s( u > V ) + J2 W ( v ) p s(u, v). (5) 
s s 

From the nature of the problem it is clear that the weight W is the same (apart from 
the change of variable). This means that it is very simple to calculate the non-nodal 
expansion. Having completed a row, we just use the updating rules in equation (J1J with 
the spreading fixed at u. This gives of a set of incomplete weights 

QsM = Y J W{u)P s {u ) v). (6) 

s 

But since Ps{u,v) is symmetric in u and v we get the desired final result as 

Pg{u,v) = Qg(u,v) + Qg(v,u). (7) 

Note that this is a very efficient algorithm. Naively one would have expected that to 
calculate the pair-connectedness one would have to sum over all 2 t realisations of the 
disorder. However, as we have just seen this is not necessary, the only complication 
being the requirement to maintain a two parameter generating function. Ultimately we 
shall calculate series for /i m>n (p), by fixing a value of a and setting u = p and v = ap. 
So to need only retain the coefficients of u l v^ provided i + j < N . The first few terms 
in the expansion for C*(u, v) are: 

CZ{u,v) = 2(u + v) 
C* 2 (u,v) = -u 4 -2u 2 v 2 -v 4 

C*(u, v) = 2u 7 - 2u 6 + (4u 5 - 6u 4 )v 2 + 2wV + (2m 3 - 6u 2 )v 4 + 4uV - 2v 6 + 2v 7 

C*(w, v) = u 12 - Au 11 + 8m 9 - 5m 8 + (2m 10 - 8m 9 + 24m 7 - 20w> 2 - (4m 8 - 8w> 3 

+ (3m 8 - 8m 7 + 24n 5 - 30m 4 )w 4 - (8u 6 - 2Au A )v b + (4m 6 - 8m 5 + 8m 3 - 20u 2 )w 6 

- (8m 4 - 24m 2 )^ 7 + (4m 3 - 4m - 5)v 8 - (8m 2 - 8)v 9 + 2mV° - 4t» n + v 12 

Note that C^(p,p) is 2*C t *(p), where C t *(p) is the non-nodal pair-connectedness of the 
pure model. In order to calculate our final series for the average cluster size S(p) and 
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Table 1. Estimates for the critical point p c and critical exponents 7, i>\\ and v± at 
various values of the disorder strength a. 
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1.0993(2) 


1.06 


0.6265485(15) 


2.2895(10) 


1.7394(3) 


1.1025(5) 


1.08 


0.6209865(25) 


2.299(2) 


1.7435(5) 


1.1068(10) 


1.10 


0.615645(4) 


2.310(5) 


1.748(2) 


1.1135(25) 


1.20 


0.59183(3) 


2.40(5) 


1.785(10) 


1.16(3) 


1.30 


0.57200(6) 


2.5(1) 


1.825(25) 


1.19(4) 


1.40 


0.55520(10) 


2.57(7) 


1.875(50) 


1.25(5) 


1.50 


0.5407(2) 


2.7(1) 


1.92(5) 


1.33(10) 


1.60 


0.5282(4) 


2.8(2) 


1.95(8) 


1.36(10) 


1.70 


0.5174(6) 


2.9(2) 


2.0(1) 


1.34(10) 



the moments /ii,o(p) and A*2,o(p) we simply fix a value of a. We then use repeated 
concatenations of C^(p, ap)/2 t to calculate the pair-connectedness C t {j>) (with a fixed), 
which in turn we use to calculate the moments. The series for the second transverse 
moment /i ,2 is obtained as /io,2 = S 2 {p) ^2 t X£(p). 

The series for C^(u, v) was derived correctly to order 115 (that is all coefficients of 
u l vi were calculated exactly for i + j < 115). This obviously results is series for S(p) 
and the other moments to order 115 as well. The algorithm used up to 5Gb of physical 
memory and the total CPU time was about 30 days. 

3. Analysis of series 

The various series were analysed using inhomogeneous differential approximants [Zj. 
Suffice to say that a .fTth-order differential approximant to a function / is a solution to 
an inhomogeneous differential equation 

K , 

J2Q t (x)(x—Yf(x) = P(x), (8) 

i=0 

where the coefficients in the polynomials Qi and P of order iVj and L, respectively, are 
chosen so that the series for the function f(x), agrees with the series coefficients of /. 
The equations are readily solved as long as the total number of unknown coefficients 
in the polynomials is smaller than the order of the series n. The possible singularities 
of the series appear as the zeros Xi of the polynomial Qk and the associated critical 
exponent \ is estimated from the indicial equation. 

Our use of differential approximants for series analysis has been detailed in 
previous papers [13 Ej and the interested reader can refer to these papers and the 
comprehensive review [Zj for further details. Typically we obtain estimates for p c and the 
critical exponents by averaging values obtained from second and third order differential 
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Figure 2. Estimates for the critical point p c and critical exponents 7, v\\ and as 
a function of the disorder strength a. 



approximants. For each order L of the inhomogeneous polynomial we average over those 
approximants which use at least the first 90% of the terms in the series. A rough error 
estimate is obtained from the spread among the approximants. Note that these error 
bounds should not be viewed as a measure of the true error as they cannot include 
possible systematic sources of error. 

In figure 121 we have plotted the estimates for the critical point p c and the critical 
exponents 7, uu and v±_ as a function of the disorder strength a (some of the estimates 
are listed explicitly in Table Q). The estimates were obtained by analysing the series 
S(p) oc - p c )~ 7 , A 4 2,o(p)/A*i,o(p) oc (p - Pc)~ Ul] and n 0j2 {p)/S{p) oc (p - p c )~ 2u± , 
respectively, and are based on results using both second and third order differential 
approximants with varying degrees of the inhomogeneous polynomial. The results are 
clearly compatible with a continuous change in all the critical parameters as the strength 
of the disorder a increases. Naturally we observe a decrease in the value of the critical 
point Pc, and in particular we note that ap c < 1. The critical exponents all increase with 
increasing a. In particular we note that uu < 2, at least while a < 1.5. It is also clear 
that the change in the values of the critical exponents is larger than the estimated error 
bars, clearly indicating that the change is not due merely to less well behaved series. 
When the disorder is relatively weak (a close to 1) the estimates for both p c and the 
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Table 2. Estimates of the critical point p c and critical exponents 7, v\\ and 2v^_ with 
a = 1.1, as obtained from third order differential approximants (L is the order of the 
inhomogeneous term). 



L 


Pc 
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v \\ 


Pc 


2v±_ 





0.6156457(10) 


2.3114(29) 


0.615644(71) 


1.83(11) 


0.6156455(42) 


2.2275(43) 


5 


0.6156454(17) 


2.3115(17) 


0.615634(35) 


1.765(23) 


0.6156474(25) 


2.2294(29) 


10 


0.6156399(16) 


2.3070(15) 


0.6156464(23) 


1.7489(11) 


0.6156449(28) 


2.2268(29) 


15 


0.6156437(22) 


2.3101(18) 


0.6156427(82) 


1.7491(30) 


0.6156458(28) 


2.2277(30) 


20 


0.6156460(26) 


2.3121(25) 


0.6156407(85) 


1.7483(26) 


0.6156447(79) 


2.2270(75) 




Figure 3. Estimates of the critical point p c and critical exponents 7, v\\ and 2j/j_ 
with a = 1.1 plotted versus the number of terms used by the third order differential 
approximants. 



exponents are still quite sharp, but as a is increased the estimated error bars become 
large and for values of a > 1.5 the exponent estimates are inaccurate. 

As a concrete example of the analysis we show some detailed results from the 
analysis with the particular value a — 1.1. In Table EJ we have listed estimates for 
p c and the critical exponents 7, v\\ and 2i/±. The estimates of p c are consistent to 5 
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digits both among the three series as well as among the approximants using different 
orders of the inhomogeneous polynomial. Likewise the estimates for the exponents are 
consistent as we vary the order of the inhomogeneous polynomial (the only exception 
being the L = and L = 5 estimates of uu). Taking into account the spread among 
the different sets of approximants (while ignoring the obviously spurious results for u\\) 
we arrive at the estimates listed in Table d To gauge whether there are pronounced 
sources of systematic errors we plot in figure H3 the estimates for the critical point p c 
and the critical exponents 7, uu and 2u± as a function of the number of terms used by 
the differential approximants. The critical point estimates are those obtained from the 
cluster size series S(p). Each point in these plots corresponds to the estimate obtained 
from a single third order differential approximant. From these plots we observe that the 
estimates do not change much as the number of terms is increased and thus conclude 
that there does not appear to be any systematic errors in the estimates. 

4. Summary 

We have introduced a simple model for directed percolation with temporal disorder, 
described a very efficient algorithm for the derivation of low-density series expansions 
for this model and presented results from the analysis of the series. The main result from 
the series analysis is that the model has a critical point p c (a) which varies continuously 
with a as does the critical exponents 7, u» and v±. The estimates for the critical 
exponent u» < 2 for most values of the disorder in apparent violation of the Harris 
criterion. 

E-mail or WWW retrieval of series 

The series for the directed percolation problem studied in this paper can be obtained via 
e-mail by sending a request to I.Jensen@ms.unimelb.edu.au or via the world wide web 
on the URL http://www.ms.unimelb.edu.au/~iwan/ by following the relevant links. 
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